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ABSTRACT 



Stellar winds are an important aspect of our understanding of the evolution of massive stars and their input into the interstellar 
medium. Here we present solutions for the velocity field and mass-loss rates for stellar outflows as well as for the case of mass 
accretion through the use of the so-called Lambert W-function. For the case of a radiation-driven wind, the velocity field is obtained 
analytically using a parameterised description for the line acceleration that only depends on radius, which we obtain from Monte- 
Carlo multi-line radiative transfer calculations. In our form of the equation of motion the critical point is the sonic point. We also 
derive an approximate analytical solution for the supersonic flow which closely resembles our exact solution. For the simultaneous 
solution of the mass-loss rate and velocity field, we describe a new iterative method. We apply our theoretical expressions and our 
iterative method to the stellar wind from a typical 05- V main sequence star, and find good agreement with empirical values. Our 
computations represent a self-consistent mass-loss calculation including the effect of multi-line scattering for an O-type star, opening 
up the possibility of applying Monte Carlo mass-loss calculations in regions of the Universe for which empirical constraints cannot 
be readily obtained. 
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1. Introduction 

Stellar winds are crucial for our understanding of the evolution 
of massive stars towards their final explosion. Furthermore, mass 
loss constitutes an important source of chemical enrichment into 
the interstellar medium. As we cannot obtain observations of in- 
dividual stars at high redshift, we need to develop a theoretical 
framework that works locally and which confidently can be ap- 
plied to the early Universe. Developing such a framework is the 
aim of this paper. 

That the winds from massive O stars are driven by radiation 
pressure through spectral lines has been known since the seminal 
papers of Lucy & Solomon (1970) and Castor, Abbott & Klein 
( 119751 hereafter CAK). Currently, there are two basic methods 
in use to compute the mass-loss rates of massive stars. These 
involve a modified version of the CAK method (Pauldrach et 
al. 119861 Kud ritzki |2002l and the Monte Carlo method (Abbott 
& Lucy 1 19831 Vink et al. 2000). The pros and cons of these two 
distinct approaches have recently been reviewed by Vink (2006). 

A common feature of supersonic flows are critical points. 
Within the CAK framework, the wind velocity, v(r), has criti- 
cal points defined by three singularity conditions, whilst assum- 
ing a power-law for the line acceleration, to obtain an approxi- 
mate solution of the equation of motion. To solve this equation, 
CAK attributed the same status to the Sobolev velocity deriva- 
tive as to the Newtonian derivative. This CAK ansatz was imme- 
di ately c riticis ed by Lucy ([19751 [19981 12007bl see also Vink et 
al. 119991 Vink 120001 . but the CAK-scaling relations have been 
widely used, presumably because of their reasonable agreement 
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with observations. The analyses of continuum emission and line 
profiles yielded mass-loss rates M. from 10~ 6 to 10~ 5 A1o/year, 
and terminal velocities of about three times the escape velocity, 
~3000 km/s for a main-sequence O-star (e.g. Howarth & Prinja 
1989). That the CAK predictions are in rough agreement with 
observations at all is largely fortuitous (see Vink 2006), as the 
assumptions of local thermodynamic equilibrium and the early 
use of limited line lists are a posteriori known to be insuffi- 
cient. The CAK framework itself however has proved to be very 
successful and improvements in non-LTE radiative transfer as 
well as atomic data have led to more reliable predictions by e.g. 
Pauldrach et al. (1986). Nevertheless, even these modified CAK 
scaling relations have not been able to reproduce the mass-loss 
rates of the denser O star winds (Lamers & Leitherer 1993 Vink 
et al. 120001 1, n or the winds of Wolf-Rayet stars (Lucy & Abbott 
[19931 Gavley [T995l Grafener et al. [20D21 Grafener & Hamann 
120051 . 

For this reason, an alternative method for predicting mass- 
loss rates, the Monte Carlo (MC) technique, was developed 
(Abbott & Lucy [T9851 Pu is [19571 Schmutz [19971 de Koter et 
al. 119971 Vink et al. 11999k The MC approach does not rely on 
the CAK ansatz, whilst still employing the computational ben- 
efits of the Sobolev approximation in radiative transfer. In the 
MC-method, mass loss does not constitute an eigenvalue of the 
differential equations that govern the stellar wind, but the rates 
follow from assuming a wind velocity structure - in the form of a 
so-called /Maw (Castor & Lamers [19791 CAK). The advantages 
of using the MC-method are manifold. There is no dependence 
on the CAK ansatz and the resulting CAK critical point, whilst 
the line acceleration is computed at all points in the wind with- 
out relying on a simplistic force multiplier approach involving 
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a 2-parameter power-law approximation for the line accelera- 
tion. Furthermore, multiple scatterings can be included naturally 
through Monte Carlo simulations. Because the kinetic energy 
density at the sonic point is negligible compared to that in the 
terminal flow, the MC-method does not rely on details of the 
trans-sonic flow (Lucy 2007a). Finally, the terminal wind veloc- 
ity is a direct and accurate (within ~10%) observable quantity in 
the local Universe. 



It is probably for these reasons that the Vink et al. (|2000) 
MC mass-loss rate£] are widely used in evolutionary calcula- 
tions, but however successful the MC approach may be in cor- 
rectly predicting the mass-loss metallicity dependence (Vink et 
al. 120011 Mokiem et al. |2007bl >, the current MC approach is 
semi-empirical. Just as we are not able to obtain empirical mass- 
loss rates from the early Universe, neither will we be able to 
obtain direct wind velocities. It is therefore important that we 
develop a reliable theoretical framework through which we are 
able to predict the mass-loss rate and wind velocity simultane- 
ously. 

We first derive analytical solutions for the general velocity 
and density structure in any mass outflow or inflow situation. 
The solutions are obtained using the so-called Lambert W func- 
tion (Muller [20011 Cranmer 2004jl. As our solutions are pre- 
sented in an explicit analytical form, physical insight can more 
easily be obtained through parameter changes. This in contrast 
to pure numerical models, where it is more challenging to pre- 
dict the response of changes in the underlying parameters. For 
the specific case of a radiation-driven wind, we do not rely on 
the CAK expression for the radiation force, rather we describe 
the line acceleration as a function of stellar radius, g{r). In ad- 
dition, the critical point of our stellar wind is the sonic point 
(as in Parker 1958 solar wind theory), and not the CAK critical 
point. The calculation of g(r) is performed through Monte Carlo 
simulations accounting for multi-line transfer, and the wind pa- 
rameters are solved simultaneously in an iterative way. 



The set-up of the paper is as follows. In Sects l2. 1142.41 the 
hydrodynamic equations for a spherically symmetric steady flow 
are introduced including a derivation of the mathematical de- 
scription of the radiative line acceleration as a function of radius 
for the case of a radiation-driven wind. The process for obtain- 
ing the fully analytical ID solutions is described and discussed 
in Sect. 12.51 Here, the velocity field for the entire family of solu- 
tions is provided in an explicit general expression from which the 
solutions for a radiation-driven wind or mass accretion flux (e.g. 
collapsing protostellar cloud) follow as unique trans-sonic solu- 
tions through the critical point. Moreover, an approximate ana- 
lytical solution for the supersonic flow is presented. In Sect. [3] 
we describe our numerical computation obtaining the radiative 
acceleration in our stellar wind models. Furthermore, a new iter- 
ative method for the determination of the consistent solution for 
the mass-loss rate is provided. In Sect.|4j we present the applica- 
tion of our models to the stellar wind from a typical 05-V-star, 
and discuss the results, before we summarise and discuss our 
findings in Sect. [5] 



2. Radiation hydrodynamics of expanding or 
collapsing systems 

2.1. Basic equations of hydrodynamics 

Considering a non-viscous, i.e. ideal fluid, the momentum equa- 
tion 



p = f - V p 

H Dt J ' 



(D 



is valid (see, e.g., Mihalas & Weibel Mihalas [T98"4"j l, where D/D t 
is the covariant Lagrangean or co-moving time derivative in the 
fluid-frame of a material element and v is its velocity, / is the 
total external body force per volume acting on a mass element 
of fluid, and V p is the divergence term of a diagonal isotropic 
stress tensor V • T, in which T = — p I and p is the hydrostatic 
pressure. 

One also needs to consider the equation of continuity 

dp 



at 



+ V - (pv) = o, 



(2) 



with the covariant divergence V ■ v of the velocity vector. 
2.2. Simplifying assumptions 

In order to solve the hydrodynamic equations analytically, and 
to account for the spherically symmetric stationary problem, we 
make the following simplifying assumptions: 

1. The stellar wind is isothermal with a temperature equal to 
the effective temperature T e ff of the central star. In this case, 
the equation of state 

p = a 2 p (3) 

is valid, where a is the isothermal speed of sound and p is 
the density of the wind. 

2. We assume a stationary ID spherical flow, 



= 0, U=f e = 



(4) 



^0. — 

dt d(p 00 

i.e. no shocks, no clumps, and no significant distortion of the 
star. 

3. In the case of a wind from a luminous early-type star, the 
wind is primarily driven by continuum plus line radiation 
forces, where the radial acceleration on a mass element is 
_ GM 

P 

with 



— a-n + «Sd 



(5) 



'rad 



(6) 



the force ratio between the radiative acceleration g™^" due 
to the continuum opacity (dominated by electron scattering) 
and the inward acceleration of gravity g, and g|'" d e (r) is the 
outward radiative acceleration due to spectral lines. M is the 
mass of the central star. 

2.3. Simplified hydrodynamic equations 
2.3.1. Wind density and mass-loss rate 

If we use the covariant derivative (see Mihalas & Weibel Mihalas 
1984) for spherical coordinates and apply assumption 2 to the 
equation of continuity we find 

d / , x . 
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for a one-dimensional, spherically symmetric and steady flow. 
Integrating this equation yields the conserved total mass flux (or 
mass-loss rate) through a spherical shell with radius r surround- 
ing the star 

M - 4nr 2 p(r)v(r) - constant , (8) 
or equivalently, for the wind density 



P(r) 



M 



4 n r 2 v (r) r 2 v (r) 



(9) 



with the defined flux F = M /4nR 2 through the star's surface at 
radius R and the dimensionless radius r = (r/R). 

Important note: All formulae derived in this Sect. [2] are ex- 
pressed in terms of f referring to the radius R, which is (through- 
out this section) the stellar (i.e. photospheric) radius of the cen- 
tral star. However, all formulae can generally also be applied 
regarding the reference radius R to be the inner boundary radius 
R m (close to the stellar photosphere), from where the computa- 
tions of a numerical wind model start (see, e.g., the results for an 
O-star in Sect.@]and associated diagrams in this Sect. 13. 

2.3.2. The equation of motion 

By using the correct contravariant components of acceleration 
(D Vj/D t) in Eq. ([TJ, for spherical coordinates, and replacing 
them by their equivalent physical components (see again, e.g., 
Mihalas & Weibel Mihalas 1984), and applying assumptions 1- 
3, we obtain the simplified r-component of the momentum equa- 
tion of a spherically symmetric steady flow 



. d 
dr 



v 



1 dp 



crit aline 

f 2 gmd p dr 



(10) 



in dimensionless form. The following dimensionless velocities 
(in units of the isothermal sound speed a) 



v 



v :- 



v, 



crit.- - \/-=-(i-r), 



a a V R 

and dimensionless line acceleration 

D 

.line . = _ .line 
orad ' ^2 6 rad " 



(ii) 



(12) 



are used. v cr jt equals the rotational break-up velocity in the case 
of a corresponding rotating star, but here it is simply the effective 
escape velocity v esc divided by a factor of V2 
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Fig. 1. The dimensionless radiative line acceleration g j.™ (r) vs. 
radial distance f (in units of R = 11 .757 Rq) in the wind from 
a typical 05-V-star (see stellar parameters in Sect. |4j. The dots 
represent the results from a numerical calculation of g 1 ™? (r,) for 
discrete radial grid points r ; . In order to determine the line ac- 
celeration parameters go, y, 6 and fo in Eq. (TBI) , these values 
were fit to this non-linear model equation resulting in the dis- 
played fitting curve (solid line) with the parameters go = 17661, 
y = 0.4758, 6 = 0.6878 and h = 1.0016 (according to 
Vk, = 3232 km/s), at the end of the iteration process (of type 
A), described in Sect. [33] and Sect. [4] respectively. 



2.4. The line acceleration term and the final equation of 
motion 

To derive a sophisticated mathematical expression for the radia- 
tive line acceleration as a function of radius r only, we collect all 
the physically motivated mathematical properties of this func- 
tion £'-(?): 

1. For hydrostatic reasons, g 1 ^ should be zero at a particular 
radius 1) near the stellar photosphere: 

= O = o, 

2. f[' a n d e should always be positive (except close to the stellar 
surface, see item 1) as the radiation from the photosphere is 
streaming outwards and so is the radiative force: 



g'Z(?)>0 for 



r > r 



Vcrit = 



V2' 



By means of Eq. (0 and applying the chain rule to the function 
1 /vt(r), we obtain 



dp 
dr 



2 1 



1 1 dv(r) 



r 3 v(f) f 2 v(r) 2 dr 
(2 1 dv(r) 

= -pin - + —— 



f v(f) dr 

Using this expression for dp/dr in Eq. ( fTOb . we finally find the 
dimensionless differential equation of motion for the radial ve- 
locity 

■ ■ 2 

(13) 



v dr 



v: * 2 

cnt , _ , ^lme 



•> rad 



that is now independent of p. 



f|.' a n d e is supposed to decrease (i.e. goes to zero) as 1/r 2 (cf. 
Castor 119741 and Abbott 1 19801 1 with increasing radial dis- 
tance f from the central star: 

Smd(r)*i>^0 for r^oo, 

g^°f has an absolute maximum somewhere in the range be- 
tween the stellar surface and the outer edge of the wind be- 
cause of the properties 1-3 above: 







All these mathematical properties for the radiative line accelera- 
tion term can be realised by the function 



^line 
Tad 



(r) 



go 

al+<5 



i-3 

r 5 



(14) 



4 



Patrick E. Miiller and Jorick S. Vink: Radiation-driven winds from massive stars 




0.99 0.995 



1 1 .005 
r / r_c 



1.01 1.015 



Fig. 2. The topology of solutions \v(r/r c )/a\ of the equation of 
motion, Eq. ( fTol i, vs. radial distance in units of the critical ra- 
dius r c , for a typical 05-V-star in the centre with the line accel- 
eration parameters go = 17661, y = 0.4758, 5 = 0.6878 and 
ro = 1.0016 in Eq. ( TBI , according to = 3232 km/s (see 
further stellar parameters in Sect. @). The horizontal line marks 
the critical velocity (i.e. sound speed v c = a). Solution 1 is the 
unique trans-sonic stellar wind solution through the critical point 
at r c = r s = 1.0110 and v(r c ) = 1.0. For the description of the 
different solutions of type 2-6, see the discussion in Sects. [2370 
and |23H 



or equivalently, 



■ 1+S (l+y) 



(15) 



which is independent of v and (d v/d f), and is a function of r 
only. 

We note that the position of the maximum (r max ,g[' a n d e (r max )) 
has to be variable and adjustable for each set of stellar and 
wind parameters, as the number and species of spectral lines 
that make up the radial force are also variable with each stel- 
lar wind. We therefore had to introduce the two parameters, go 
and y. Moreover, item 1 is fulfilled at radius f' = fJ s and item 3 
can be guaranteed by appropriate values of the exponents 5 and 
y (where < y < 1 and < 6 < 1). Altogether, we had to intro- 
duce at least the set of four parameters in Eqs. (TBI and (15[ . 

Thus, the equation of motion ([TBI becomes 



r 



dr 



v = — 



cut 



So 



2 



(16) 



This equation is fully solvable analytically as we show later. 

2.5. Analytical solutions of the equation of motion 

2.5.1. The critical point and critical solutions 

The equation of motion (Eqs. [T3l and [TBI) yields several families 
of solutions that have quite different mathematical behaviour and 
physical significance (cf. Fig. [2j. 

The left hand side of Eq. ( fT6b vanishes for [dv/dr + Q)f c at 
the critical radius r c , where 



v (fc) = v (f s ) = 1 



(17) 



That is, the critical point velocity here is equal to the isothermal 
sound speed v = 1, and the critical radius is just the sonic radius 



as is also the case for thermal winds or mass accretion events 
(where £™ * 0). 

We are now interested in under what conditions one can ob- 
tain a continuous and smooth trans-sonic flow through the criti- 
cal point r c of Eq. ( fT6] >. For the case of a stellar wind, this means 
how to obtain a smooth transition from subsonic and subcriti- 
cal flow (v < v c = 1) at small r < r c to supercritical and super- 
sonic flow (v > v c ) at large f > r c , when this critical solution has 
a finite positive slope (dv/dr) > at f — f c (cf. solid curve 1 in 
Fig. IJji'O Then, it is evident from the left hand side of Eq. ( fTST ) 
that one can obtain such a trans-sonic wind, if the right hand side 
(1) vanishes at the critical radius r c , (2) is negative for r < r c , and 
(3) is positive for r > r c . 

The opposite situation occurs for the case of mass accretion 
in e.g. a collapsing cloud. If (dv/dr)? c < 0, we obtain the second 
unique trans-sonic and critical solution in which v (r) is mono- 
tonically decreasing from supersonic speeds for r < f s , e.g. near 
the protostar, to subsonic speeds for r > r s at the outer edge of 
the cloud (see also the second solid line 2, in Fig. [2] for the case 
of a corresponding accretion flow with a star as the central ob- 
ject). 

Here we are mainly interested in the critical wind solution of 
Eq. JTol l. The right hand side of Eq. ( [Tol l vanishes at the critical 
radius f c that solves the equation 



2r^^>-e it ^ (1+r) - I +fo(^-r ) y =0. 



(18) 



Therefore, the critical radius has to be determined numerically 
by means of the above equation and the line term parameters go, 
y, 6 and ro- However, if one assumes values of y and 6 close to 
1, one can provide a good approximation (i.e. analytical solution 
of Eq. (fT8l ) for the critical radius 



(19) 



For the special case of a thermal wind, where go can be set to 
zero, we obtain 



cnt 

T 



GM(\ -D 
2aR 



(20) 



In some special cases of mass accretion (e.g. a collapsing cloud), 
one can even set F = 0. 



2.5.2. Solving the equation of motion 

The Equation of motion ( fT6] i can be solved by first integrating 
the lefthand side over v, and then integrating the righthand side 
over r, separately, which yields 



v - In v = 2 + 4 lnr + 

f r 8 (l+y) 



l+y 



l-£ +C, (21) 



with the right hand side of Eq. ( l2Tl i denoted as the function 



f{r, f,v):= 2^+41nr+ 2 
r 



ro 



. . ,1-^ +C(r', v')(22) 

ro6(l+y)\ r° ' 



l+y 



with the constant of integration C that is determined by the 
boundary condition of the radial velocity v' at a given radius 
f'. 



2 We furthermore assume that both v and (dv/dr) are everywhere 
single- valued and continuous. 
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2.5.4. The solution(s) of the equation of motion 

It is now possible to provide an explicit analytical expression 
for the solution v of the equation of motion ( fToT i. or Eq. (l25t . 
by means of the W function. If we compare Eq. d25l l with the 
defining equation of the Lambert W function (Eq.l26l>. 

Wi(jc) e WlW = x , 

we find that 

-v 2 = W k (x) 

or 

v = ± V-W t (x) (27) 



Fig. 3. The two real branches of the Lambert W function: Wo(i) 
(full curve) and W_i(x) (broken curve). 



is the general solution of the equation of motion that passes 
through the point (f' , v'), with the argument function of the W 
function 



x(r) 



a -/(r;r',0') 



(28) 



From Eq. (EB, we can determine C (?' , V) for the solution Since the argument of the W function in Eq. (E2J is always real 

and negative, it is guaranteed that the argument of the square 
i +r root never becomes negative, and hence the solution is always 
.(23)"eal (see Fig. [3j. 

Inserting / (f; ¥ , v') from Eq. d24l > into Eq. 



that passes through the particular point (r', v'), 

i\2 



C(r',v') = i>' 2 -lnv' 2 -2 



V • 2 

cnt A , ~/ 

4 In r - 



go 



r ?q 6 ( 1 + y) 

Therefore the function / in Eq. ( l22l becomes 

/ (r; f', v') = v' 2 - In v' 2 + 2 v> 2 nt ij - i) + 4 In [L 



ro_ 



yields 



2 

+ — 



go 



i-3 

r 5 



i+r 



(24) 



fo 6 (1 + y) 
And from this, Eq. (f2TT > now reads 
v 2 _lnv 2 = /(r; r,v') 
or, equivalently, 

- v 2 e- p2 = - & -fW' y \ (25) 

which is solved explicitly and fully analytically in terms of the 
Lambert W function. 



x(f;f', v') = - I — ) v exp 



So 



r 5(l +y) 



-2v 



!-3 



(i 


1 


) 




~ F 




i+ r 








('- 









l+y\ 



the general expression for the argument function x depending on 
the parameters (f' , V). 

Thus, for the trans-sonic case of a stellar wind or general 
accretion flow, where f — f c = r s and V — v c = 1, the analytical 
solution is 



v(r) = ±fW)) 
with 



(30) 



2.5.3. The Lambert W function 

The W function is defined to be the implicit function satisfying 

W(z)e wfe) =z (26) 

(cf. Corless et al. FT993lfT996b . 

As the equation y(x)expy(x) = x has an infinite number of 
solutions y(x) for each non-zero value of x, W has an infinite 
number of branches. We are only interested in the physically 
relevant case, where x is real and -1/e < x < 0. Then, there are 
two possible real values of W(x) (see Fig.0. As one can already 
see from the defining equation of W, the W function vanishes 
at x = 0, is negative for x < 0, positive for positive values of 
x, and must be -1 at the point x = -1/e. The branch satisfying 
-1 < W(x) for x in the range of [-1/e, oo) is denoted by Wq(jc) or 
just W(x), and the branch satisfying W(x) < — 1 over the inter- 
val [-l/e,0) by W_i(x). The branches Wo(x) and W_i(x) share 
the branch point at x = -1/e. Wo(x) is referred to as the prin- 
cipal branch of the W function, which is the only branch that 
is analytic at 0. The other remaining non-principal branches of 
W all have a branch point at 0, and they are denoted by Wfc(x), 
where k is a non-zero integer. 



x(f) = - 



exp 

go 



-2v 



f S(l +y) 



1-3 

f 5 



1 1 

f f t 

l+y 



- 1-^ 



l+y\ 



- 1 



(3D 



From this equation, one can easily derive the solution for the 
simpler cases of a thermal wind (e.g. the solar wind), or for a 
particular accretion flow, as e.g. a collapsing molecular cloud, 
where go can be set to zero and v 2 = 2 r c (cf. Eq. l20t: 



v(r) = 



-w* 



(32) 



Here, the minus sign (in front of the square root) refers to the 
collapsing system. 

We are only interested in the possible two real values of 
W(x), the k = 0, -1-branches in Eq. d27b or Eq. d30l l. where x 
is real and — 1/e < x < 0. The branch point at x = -1/e, where 
these two branches meet, corresponds to the critical point f c , 
where the velocity in Eq. (130b becomes v = v(r c )= 1 . Depending 
on which branch of W one is approaching this point x = -1/e, 
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Fig. 4. The argument function x(r) of the W function, Eq. Pit , in 
the trans-sonic velocity law, Eq. O0I) . vs. radial distance f = r/R 
for the stellar wind from a typical 05-V-star with the line ac- 
celeration parameters go = 17661, y - 0.4758, 6 = 0.6878 and 
ro = 1.0016 in Eq. (1141 1. according to Vco = 3232 km/s (see fur- 
ther stellar parameters in Sect.|4j>. The minimum of this function 
is at the critical radius r c = r s = 1.01 10 where x(r c ) = -1/e. 



one obtains a different shape of the v(r)-curve, i.e. a stellar wind 
or a collapsing system. 

However, to determine which of the two branches to choose 
at a certain range of radius between [l,r c ] and [r c ,oo) so as to 
guarantee a continuous, monotonically increasing, and smooth 
trans-sonic flow (as in the case of a wind), one still needs to 
examine the behaviour of the argument function x(r) of W, in 
Eq. ( l3T1 l. with radius f (see Fig. |4j». 

The argument function x(f) decreases monotonically from 
the stellar radius f — 1 (with a value of nearly zero) to its mini- 
mum at f — f c with x = -1/e to afterwards increase monotoni- 
cally again. 

Stellar wind: By choosing first the real principal branch 
Wo(x) of the Lambert function that is followed in the negative 
direction until the critical point x(r c ) = -1/e and Wo(x)= -1 
(i.e. v(r c )) for rising radii from near the star surface f « 1, and 
hence coming from x(r) < and Wo(x)< 0, yields the desired 
critical wind solution. This choice causes the strict monotonic 
increase of the velocity v(f) from nearly zero until the sonic 
point (in the subsonic region). To ensure that the solution passes 
smoothly through this critical point, i.e. that the derivative v'(r c ) 
is continuous, one has to change there to the W_i(x)-branch. 
From there onwards, the real part of this branch is followed fur- 
ther in the negative direction W_i(x)— > - oo for larger radii and 
therefore with decreasing \x(r)\ — > 0, again. From this, one ob- 
tains the desired continued monotonic slope of v(f) in the super- 
sonic region for f > r c (= r s ). 

Accretion solution: In this case, the unique critical and 
trans-sonic solution \v(r)\ is strictly monotonically increasing 
from small values in the subsonic region far away from the cen- 
tral object (for f > f s and |v(r)| <K 1), into the supersonic re- 
gion (for f < r s and \v(r)\ > 1) with declining radius f. Here the 
slope v' is negative and continuous at the sonic (critical) radius r s 
(as a necessary boundary condition). This case can be achieved 
when we again first take the Wo(x)-branch for > x(r) > —1/e 
(in the subsonic region) and then change over at the sonic point 
to the W_i(x)-branch for -1/e < x(r) < 0; now however with 
decreasing radius f — > 1, Here one is again wandering in the 
same negative direction through the possible two real values of 
W(x) in the range of - 1 /e < x < 0. We are approaching the crit- 



ical point x(r s ) = — 1/e (where the argument function x(r) has its 
minimum, cf. Fig. |4]i with declining radii r > r s from the right 
side and from x(f) < 0, in order to change the branch at this 
point, where x(r) afterwards increases again strictly monotoni- 
cally (with smaller radii) until almost x(r + ) < near the central 
object. 

To summarise, the velocity in the radial direction (i.e. the 
exact analytical and spherically trans-sonic solution of our equa- 
tion of motion) can be described, 

(a) for the case of a stellar wind, by Eq. (1301 (and the posi- 
tive sign in front of the root) and the argument function in 
Eq. (l3T1 >, choosing the branch 



k = 



for 1 < f < r c 
-1 for f > r c 



(33) 



of the W function at a certain radius f, and 
(b) in case of a general accretion flux, as well, by Eq. (l30l (but 
now with the negative sign in front of the root) and the argu- 
ment function in Eq. dBTt , choosing the branch 



-1 for 1 < f < r c 
for r > h 



(34) 



(O 



depending on the location r, and 

in the special cases of a thermal wind or collapsing system 
like a collapsing protostellar cloud, by Eq. ( |32l , choosing 
the same branches as mentioned above in case (a), or (b) 
respectively, and the appropriate sign (in front of the root). 

In addition to these two critical solutions (type 1 and 2, cf. num- 
bering in Fig.|2]l, already discussed, that pass through the critical 
point (i.e. sonic point), all the other four types of solutions were 
obtained from our general velocity law, Eq. ( l27l i with Eq. (|29] >, 
choosing the following branches of the W function and values of 
(f', V), for the point we demand the solution to go through: 

- Type-3: Subsonic (subcritical) solutions 
k = , f = r c , v' < 1 

- Type-4: Supersonic (supercritical) solutions 
k = -1 , f = f c , v > 1 

- Type-5 & 6: Double-valued solutions 
k — and k — -1 for 

f' , v' = arbitrarily, where f' +r c , v' + 1 . 

Type-3 solutions are everywhere subsonic (choosing only the 
principal branch, k — 0, of the W function). Those of type 4 
are everywhere supersonic (selecting only the k = -1 -branch in 
the velocity law), and those of type 5 and 6 are double-valued, 
composed of both the k — and k — - 1 -branch, below and above 
the sonic line, respectively. In this connection, the two sub- and 
supersonic pairs of curves of this last mentioned types, subdi- 
vided into (5a, 6a) and (5b, 6b) in Fig. [2] belong together. They 
are fixed, not only by the same chosen branch of W in Eq. ( |27| >. 
but also by the same selected parameters for the solution through 
the identical given point (r' , V). 

From our analytical wind solution, we find that v increases 
without limit as f — > oo. This unphysical increase in v at large f 
is an artificial consequence of our assumption (1) that the stellar 
wind is isothermal with a temperature, T e ff , of the central stellar 
photosphere. This assumption is valid close to the photosphere, 
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Fig. 5. Comparison of the approximated wind solution (invisible 
dashed line, overlaid by the solid line) with the exact solution 
v (r) (in units of sound speed) vs. radial distance f (in units of R = 
1 1 .757 Rq) for the wind from a typical 05-V-star with the line 
acceleration parameters go = 17661, y = 0.4758, 6 = 0.6878 
and ro = 1.0016 in Eq. ( fT4l >, according to Voo = 3232 km/s (see 
stellar parameters in Sect. 0}. The approximated curve was ob- 
tained from Eq. d39l l, where the exact solution was obtained from 
Eq. d30l . together with Eq. PTT l and Eq. (1331 . with a critical ra- 
dius of r c = 1 .01 10. For a more detailed comparison see Fig. Q 
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Fig. 6. Comparison of the approximated wind solution (dashed 
line) with the f3 wind velocity law (solid line) v(r) (in units 
of sound speed a) over a large radial distance f (in units of 
R — 1 1 .757 Rq) for the wind from a typical 05-V-star. The ap- 
proximated wind curve was obtained with the same parameters 
as in Fig. [5] where the /3 wind velocity curve was retrieved from 
Eq. (|4TT i with an exponent of (3 — 0.7379, a terminal velocity of 
(voo/a) = 177.9 (according to v M = 3232 km/s) and a parameter 
value of r' Q = 1.0095 (used in IS A- Wind). Both curves approach 
the same velocity limit Voo as f goes to infinity. 



but may not hold at large distances, where the gas expansion will 
ultimately force the wind to cool down. By insisting on keep- 
ing the wind temperature T fixed, we introduce an undesirable 
energy source to do work on the gas and accelerate it without 
limit. In the case of the solar wind, this classical problem was 
overcome in the work of Parker (1958) (see also the review of 
Mihalas 1978) by assuming the corona near the sun to be isother- 
mal, yet allowing the wind to expand adiabatically at larger radii. 

We overcome this problem by deriving an analytical expres- 
sion for the wind solution in the supersonic approximation, in 
the following Subsect. 12.5.51 by neglecting the pressure term in 
the equation of motion which introduces a terminal velocity v M , 
as well as a relationship between v„ and our wind (and line ac- 
celeration) parameters go, y, 5 and ro. 

At the end of the iteration process (described in Sect. l3.3t . we 
obtain the converged values of the previously mentioned wind 
model parameters, together with the real position of the criti- 
cal radius r , to be able to evaluate the final wind solution with 
Eqs. ( 1301 ) and d3TT >. This way, the complete (exact) wind velocity 
curve is produced (in Fig. [5]) to be compared with the approx- 
imated wind solution over a wide wind range from the stellar 
photosphere up to a distance of 20 stellar radii. This compari- 
son (see also the enlarged images of Fig. [7]) shows a high level 
of agreement, and only a weak discrepancy between both curves 
v (r), in particular for larger radii, where the (even) unlimited 
function for the exact wind solution is close to the 'speed' lim- 
ited approximated solution curve. Through this way of adjust- 
ment, the undesirable unlimited mathematical behaviour of our 
expression for the wind velocity has become insignificant over 
the whole r range of interest. 

The subsequently derived approximated wind solution is 
only valid in the supersonic region and is not supposed to be ap- 
plied to the subsonic region, where it even becomes imaginary, 
particularly in our wind model in the range of < r < r s , see 
upper diagram of Fig. [7] This explains the only noticeable dis- 
crepancy between both wind velocity curves in that radial range, 
from close above the sonic radius down to the stellar surface. 



The numerically computed model atmospheres that we ap- 
ply later in our iteration process allow a temperature stratifica- 
tion T (r) of the wind, i.e. we do not assume a fixed T as in our 
analytical expressions. 

2.5.5. Approximated solution of the equation of motion 
By neglecting the pressure term in the Eq. of motion ( [Tot . 

1^0, 
p dr 

which is a good approximation for the stellar wind solution in 
the supersonic region with f > r c = r s , Eq. ( fTBI l becomes 



dr 



i 9 = _!£rit + ^0 

a2 sl+<5(l+y) 



=5C-*)' 



(35) 



This simplified equation of motion can again be solved by first 
integrating the left hand side over v, and then integrating the right 
hand side over f, separately, which yields 



~2 V 'crit ^ 2 g0 

v = 2 1 

r r 6 (1 + y) 



i-5 

fS 



1+)' 



C. 



(36) 



In order to determine the integration constant C, we assume a 
boundary condition v (?) « for the wind velocity at radius f = 
r l J 5 , close to the stellar photosphere, i.e. 



C{r =r ' ,v 



= 0) 



v 2 . 

cut 



(37) 



Thus, Eq. d36l ) now reads 



v 2 = 2 



1 



M(i +r) 



l+y\ 



(38) 





17.5 




15 




12.5 


CO 




£ 


10 






7.5 


> 






5 




2.5 





147.5 




145 




142.5 


CO 






140 


£ 




i 


137.5 


> 






135 




132.5 




130 



160 

159 

* 158 
£ 

| 157 
156 



Patrick E. Miiller and Jorick S. Vink: Radiation-driven winds from massive stars 

2.5.6. Comparison with the ft velocity law 
Eq. ( l39b yields a terminal velocity v„ (as r — > oo) of 
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Fig. 7. Model results for the wind velocity v(r) (in units of 
sound speed a) vs. radial distance f (in units of R — 1 1.757 Rq) 
from a typical 05-V-star with the line acceleration parameters 
go = 17661, y = 0.4758, 6 = 0.6878 and r = 1.0016 (cf. 
Eq. IT4b . and a critical radius of f c - f s - 1.0110 (see stellar 
parameters in Sect. @). The approximated (dashed line) and the 
exact wind solution (solid line) are compared with each other 
for three different and magnified radial ranges: Very close to the 
central star between r — [1.0, 1.1] (see upper diagram), at inter- 
mediate distance of f = [5.0, 10.0] (see middle diagram) and far 
away between f = [15.0, 20.0] (see lower diagram). 



Vco = -I — 




S(l+y) 



_~2 fl-l/S 
V crit r 



(40) 



which is now comparable to the Voo parameter in a (so-called) j3 
velocity law (cf. Castor & Lamers 1979, CAK) 



v (r) = Voo 1 1 - — 



(41) 



To be able to compare the y (and 5) parameter in our wind law 
with the exponent in the f3 law (as we use ft as an input param- 
eter in our model atmosphere calculations), we express our line 
acceleration parameter go in terms of Voo by means of Eq. d40l , 
i.e. 



I "2 
V 



V 

— + 

2 r. 



2 > \ 
crit 



'0 ' 

and insert it into Eq 



r (5(l +y) , 

as to obtain 



(42) 




si-i/a 



tl. 2 

2 V ~ 



cnt 



l+y\ 



which now depends on v' m , y, 5 and ro. 

Since 5 is of the order of 1, as is ro, one can approximate 
the expression ri -1 ^ 5 in Eq. ( l43l l as 1. Furthermore, as our line 
parameter ro is defined as the parameter for which the line accel- 
eration becomes zero at radius r = r l J s (cf. Eq. fT4b . this radius 
is very close to the radius r' (j in the /Maw (in Eq. HTt . where the 
wind velocity is assumed to be zero, i.e. * ro. 

Then, we can set the velocity in Eq. ( f4TT > equal to our ve- 
locity law in Eq. d43l >. to search for a relationship between the 
parameters f3, y and 6, which yields 



1-2* 



2/3-1 



-b + (\+b) 



!-3 



(44) 



with 

ro \ Vc, 



Herein, for large radii f (and especially for small values of b, 
e.g. b < 0.1 for an O-V-star), the right hand side can be approx- 
imated by the last fraction only, which leads to 



or equivalently 



(45) 



resulting in the approximated wind solution 




v% I — 

cnt \ J. 



go L _ [o 
6(l+y)\ f 6 



which can be expressed without the W function. 



l+y\ 



(39) 



2/3 
1+7 



log 1 - 



(46) 



Next, the righthand side in Eq. ( l46b can be approximately set to 
1, for values of 5 — > 1 or generally for smaller distances from 
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the central star in the supersonic region as f - . , (| 
results the desired relationship between ji and y 

2j3~ 1 +y, 



IfS 



1. This 2.5.9. The reasons for deriving the approximated solution 



(47) 



independent of 6, that is most valid for the previously mentioned 
values for 6 at smaller radii r. It also applies at very large dis- 
tances f » 1, since then, the numerical values inside the brack- 
ets of Eq. ( |45t are close to 1 and this equation is fulfilled for any 
value of the exponents ft and y in all cases. Only for intermediate 
distances from the star at lower values of 6 (not close to 1), the 
relationship between /? and y is possibly not well approximated 
by Eq. WT\ . 

Figure [6] illustrates this described behaviour by the compar- 
ison of the approximated wind solution with the /3 law for our 
wind model calculations from a typical 05-V-star. 

2.5.7. Fitting formula for the line acceleration 

Thus, we can provide another expression for the line acceleration 
(in Eq. [T4"t . now dependent on and on the parameters /3 (or y 
respectively), 6 and ro (by eliminating go using Eq.l42b: 





( "2 

2 

V 


. cnt 

r o > 


or 






&) = 
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V 


2v 2 ,) 

, cnt 
r > 



r 6(l +y) 



1 

t-1+6 



1 - 



l 



!-8 



2/3-1 



(48) 



(49) 



Both of these non-linear expressions can then be used as fitting 
formulae and applied to the results from a numerical calculation 
of g|'" d e (r,) for discrete radial grid points r,, in order to determine 
the line acceleration parameters y (or equivalently /3), 6 and ro, 
and the terminal velocity X (cf. Fig. [TJ. 



2.5.8. Physical interpretation of the equation of critical radius 

Through the use of the exact wind solution, by using Eq. ( fT8l . 
valid for the critical radius r c , we can solve for the line acceler- 
ation parameter go and insert it into Eq. d40b from the approxi- 
mated wind solution, to provide another expression for the termi- 
nal velocity depending on the location of the critical (i.e. sonic) 
point 



v M = 



fi-h 



51(5-1 



6(1+7) 



fet- 2 ?c)- 



-2 -1-1/5 

v ,r„ 

cnt u 



(50) 



Or vice versa, by Eq. ( fT8l , the location of the critical point 
(through which the exact analytical wind solution of our equa- 
tion of motion EOM ([Tol l passes) is mainly determined, on the 
one hand, by the given terminal velocity v M , via the line accel- 
eration parameter go- On the other hand, the position of r c must 
also be dependent on the given minimum velocity v m at the inner 
boundary radius R m , where the velocity solution passes. This in- 
ner velocity v; n follows indirectly from the other remaining line 
acceleration or wind parameters y, 6 (which make up the shape 
of the velocity curve) and especially ro (where the value of the 
latter parameter is determined by the radius r^ 6 at which g'™ e is 
zero). 

Since the inner boundary condition of the velocity v; n is con- 
nected to the mass-loss rate M, through the equation of mass 
continuity by Eq. ([8]) and the given density at the inner bound- 
ary, the position of the critical radius is uniquely specified by the 
values of v' m and M. 



The reason for the derivation of our approximated wind solution 
(in Sect. 12.5.51 ). that is only valid in the supersonic approxima- 
tion, is not its application to represent the unique wind solution 
in the whole wind regime (i.e. including the subsonic region). 
The approximated solution is required to be able to calculate the 
terminal velocity belonging to our exact analytical wind so- 
lution described by Eqs. fl30l ) and d3T1 > and to overcome its arti- 
ficial and unlimited increase with increasing radius by assuming 
an isothermal wind (as discussed at the end of Sec. |2.5.4T >. 

In other words: the approximated solution is necessary to be 
able (1) to find the relation between Voo and the line acceleration 
(or wind) parameters go, y, 6 and ro, and (2) to compare the 
exponents y and 6 to the exponent/? (in the /J-law). 

This, in turn, is necessary because our applied model atmo- 
sphere code assumes a /3-velocity law for the whole supersonic 
wind region which is mainly determined by the input parameters 
Vco and /3. The additionally used code parameter r' (cf. Eq. |4TT ) 
is related to the velocity v m at the inner boundary radius which, 
again, is directely dependent on the prescribed mass-loss rate M 
via the continuity Eq. (O (and the density at the inner boundary 
given therein). 



3. Numerical methods 

3.1. Computing the radiative acceleration 

The calculation of the radiative acceleration of a stellar wind re- 
quires the numerical computation of the contribution of a very 
large number of spectral lines. We first calculate the thermal, 
density and ionisation structure of the wind model by means 
of the non-LTE expanding atmosphere (improved Sobolev ap- 
proximation) code ISA-Wind (de Koter et al. 1 1993L [T997l ). As a 
next step, we calculate the radiative acceleration as a function of 
distance by means of a Monte Carlo (MC) code MC-Wind (de 
Koter et al. I1997I Vink et al. |1999) . The basic ideas behind this 
technique were first applied to the study of stellar winds from 
early-type stars by Abbott & Lucy dl985l l, but in our calculation 
we account for the possibility that the photons can be scattered 
or absorbed and re-emitted (due to real absorption) or eliminated 
(if they are scattered back into the star). 

The radiative transfer in MC-Wind involves multiple contin- 
uum and line processes using the Sobolev approximation. The 
continuum processes include electron scattering and thermal ab- 
sorption and emission, whilst the line processes include photon 
scattering and photon destruction by collisional de-excitation. 
Line processes can only occur at specific points in each shell 
of the stellar wind, whereas continuum processes can occur at 
any point. To decide whether a continuum or line process takes 
place, we apply a similar approach to Mazzali & Lucy ( 11993I I, 
where a random optical depth value is compared to the combined 
optical depth for line and continuum processes along a photon's 
path. Our code was improved in the following way: after it is de- 
cided that the next process will be a continuum process, a second 
random number is generated to decide which continuum process 
takes place, an electron scattering or absorption process (as im- 
plemented in Vink et al. [19991 Vink|2000]l. 

The radiative acceleration of the wind is calculated by fol- 
lowing the fate of the large number of photons where the atmo- 
sphere is divided into a large number of concentric thin shells 
with radius r and thickness dr, and the loss of photon energy, 
due to all scatterings that occur within each shell, is determined. 
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The total line acceleration per shell summed over all line scat- 
terings in that shell equals (Abbott & Lucy 1985) 



M dr 



(51) 



where -dL (r) is the rate at which the radiation field loses energy 
by the transfer of momentum of the photons to the ions of the 
wind per time interval. 

The line list that is used for the MC calculations consists 
of over 10 5 of the strongest lines of the elements from H to 
Zn from a line list constructed by Kurucz (1988). Lines in the 
wavelength region between 50 and 10000 A are included with 
ionisation stages up to stage VII. The number of photon packets 
distributed over the spectrum in our wind model, followed from 
the lower boundary of the atmosphere, is 2 x 10 7 . The wind is 
divided into 90 spherical shells with a large number of narrow 
shells in the subsonic region and wider shells in the supersonic 
range. 



3.2. Computing the mass-loss rate for known stellar and 
wind parameters 

By neglecting the pressure term and using the expression for the 
line acceleration per shell (Eq. |5TV an integration of the Eq. of 
motion (fTDl >. from stellar radius to infinity, yields (cf. Abbott & 
LucyQ985) 




^(vi+vL) = AL, 



Fig. 8. The numerically computed wind velocity (assumed by 
ISA- Wind, see solid line) in comparison to the exact wind solu- 
tion (dashed line) and the /3 wind velocity curve (dotted-dashed 
line), also used in ISA-Wind, vs. radial distance r (in units of 
R = 11.151 Rq) close to the 05-V central star (in Sect.!). The 
exact solution fulfills the Eq. of motion ( [ToT l for the converged 
parameter values go = 17661, y = 0.4758, 5 = 0.6878 and 
ro = 1.0016 (shown in the last row of Table|2|, at the end of the 
iteration process described in Sect. 13.31 In this case of conver- 
gence, the critical radius r c of the Eq. of motion ( fTSI ), determined 
by Eq. ( fT8l ). becomes equal to the sonic radius r s = 1.01 1 deter- 
mined by ISA- Wind (as demanded by our theory and iteration 
method). At this radius f = 1.011, both wind velocity curves, 
i.e. the exact and the computed velocity curve by ISA- Wind, in- 
tersect each other. 



or equivalently, 



M = 



2AL 



vt, + K 



(52) 



where AL is the total amount of radiative energy extracted per 
second, summed over all the shells. This equation is now fun- 
damental for determining mass-loss rates numerically from the 
total removed radiative luminosity, for the prespecified stellar 
and wind parameters v esc and v„, respectively. 

3.3. The iteration method: Determination of M and the wind 
parameters 

In computing the mass-loss rate and the wind model parameters 
from a given central star with fixed stellar parameters L, T e ^, R, 
M, T, the following iterative procedure is applied: 

1 . By keeping the stellar and wind parameters M„, Voo n , [}„ vari- 
able throughout our iteration process, we use arbitrary (but 
reasonable) starting values M_i, v m i , f3-\ in iteration step 
n = -1 (cf. Tables E-0. 

2. For these input parameters, a model atmosphere is calculated 
with ISA- Wind. The code yields the thermal structure, the 
ionisation and excitation structure, and the population of en- 
ergy levels of all relevant ions. Then, the radiative acceler- 
ation gj.™ (ji) is calculated for discrete radial grid points r,- 
with MC-Wind and Eq. (1511 . In addition, an improved es- 
timate for the mass-loss rate M° ut is obtained by Eq. d52b . 
which can be used as a new input value for the next iteration 
step. Moreover, one obtains a new output value for the sonic 
radius r Sn (which has to be equal to the critical radius r c of 
our wind theory). 



To determine the improved line acceleration parameters y„ 
(or equivalently /?„), 5 n and fo n , Eq. (|48l or d49l is used as the 
fitting formula to apply to the numerical results for g'™ (r,), 
cf.Fig.CD 

By applying Eq. < f30b and inserting the current values of pa- 
rameters y,„ 6„ and ro„, as well as the current sonic radius r Sn 
for f c , we obtain a new approximation of the terminal veloc- 
ity Voo n , i.e. 



<J n (i +r«) 



( v 2 . — 2 r s ) - v 2 .,r\ 1 

V cnt s n / cnt ()„ 



5. With these improved estimates of M n , Voo„, Pn as new input 
parameters, the whole iteration step, defined by items 2-4, 
is repeated until convergence is achieved. 

3.4. Physical explanation of the iteration method for 
calculating the unique wind solution 

In ISA- Wind the equation of motion is not solved explicitly for 
the supersonic wind region and, additionally, under the assump- 
tion of any input values of the stellar and wind parameters M, Voo 
and y6, any arbitrary 'solution' for the wind velocity field is con- 
ceivable. However, as from the theoretical point of view, it must 
be r c = r s . Hence, to find numerically the unique wind solution 
by application of ISA- Wind, Eq. (l53l l guarantees that at the end 
of the iteration process this condition f c = r s is just fulfilled for 
the sonic radius r s in ISA- Wind for the converged set of wind 
parameters. 

Eq. d53l l follows from this condition for the critical radius 
inserted into Eq. ( fl~8b . Since the line acceleration parameter go 
can be expressed by means of the remaining parameters (y, 6, ro) 
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via the terminal velocity v M (cf. Sects. |2"378l and |2~5.9l l, Eq. d40l > 
has additionally been engaged for its derivation. 

Thus, Eq. ( l53l determines the terminal velocity v Mn that our 
exact wind solution of the EOM ( fTol ) would need to have, in 
order to pass through the critical point at radius r Cn = r Sn for 
the parameters of the fit to the numerical line acceleration curve 
obtained previously in the same iteration step (y„, 6„, ro„). This 
improved estimate of Voo„ is then used as a new input parameter 
in ISA- Wind for the terminal velocity in the next iteration step. 

Naturally, the (current) position of the critical point is also 
affected by the parameter value of M„, since the value of M is 
simultaneously being determined iteratively and numerically to- 
gether with the other remaining wind parameters by MC-Wind 
and ISA- Wind (with its assumed velocity- and density field) in 
our iteration process. By the improved estimate of M at each it- 
eration step and the use of ISA- Wind, the lower boundary wind 
velocity V[ n (at the inner boundary radius) is also gradually being 
adapted to an improved value (via the continuity Eq. ([8j and the 
given inner boundary density). 

Finally, at the end of the iteration process in the case of the 
convergence of all wind parameters, both different wind velocity 
curves (i.e. the exact solution of EOM ( fTol i and the wind curve 
computed by IS A- Wind) share the same sonic point, i.e. intersect 
each other at radius r s (see Fig. |8j, in order to approach the same 
terminal velocity value Vco at large distances from the central star 
as r — > oo (cf. Fig.[6]l. 

Then, our exact wind solution, analytically expressed by 
Eqs. d30j and CO), fulfills the EOM (O for the converged wind 
parameters go (or v«,), y, 5, ro. It now represents the unique wind 
solution through the critical point of the EOM ( fTol i for the wind 
from a given star with fixed stellar parameters L, T^, R, M, T. 

3.5. The model assumptions and velocity field in the 
atmosphere code ISA-Wind 

The ISA- Wind code has been described in detail by de Koter et 
al. ([1993 1996). Here, we will only discuss those model assump- 
tions in IS A- Wind which affect our wind formalism by using the 
code in our iteration process. 

In ISA- Wind the wind velocity structure is divided into two 
regions: 

1. The low-velocity subsonic region determined by an EOM, 
similar to our Eq. ( fTol i. however in contrast to this, neglect- 
ing the line acceleration term, and, with an additional term 
oc dT/dr on the rhs of the Eq. (assuming a non-isothermal 
wind), and 

2. The high- velocity supersonic region where the ft velocity law 
in Eq. (l4TT > is assumed to represent the wind solution of the 
EOM (i.e. including then the line radiation forces). 

These two regions are connected at a radius f con , which is deter- 
mined, together with the parameter r' Q (in Eq. |4TT >. with the con- 
dition that the velocity law is continuous in v and dv/dr at this 
connection radius. The value of r con is generally located close to 
and below the sonic point r s (cf., e.g., Fig. [8] where r con = 1.010, 
r = 1.0095 and r s = 1.011). 

The inner boundary condition of the flow velocity, Vi„ at the 
inner boundary radius R m , follows from prescribed (i.e. freely 
selectable) values of the density at the inner boundary condi- 
tion and mass-loss rate (through the equation of mass continuity 
in form of Eq. [8j. Typically Ri n is chosen to be situated at a 
Rosseland optical depth of tr > 20. Therefore, Ri n does not cor- 
respond to the stellar (i.e. phostospheric) radius R, that depends 



significantly on frequency. In IS A- Wind the photospheric radius 
Rphot is defined as the point where the thermal optical depth at 
5555A equals r e = 1/ V3. In practice is R m < /? p hot- 

Hence, ^ p h t follows from the model computation and the 
chosen model parameters describing the density and velocity 
structure in the photosphere: the radius R m and density p^ at 
the inner boundary. The wind is specified by the input param- 
eters M, Vco and /3. The parameters specifying the temperature 
structure are the luminosity L, R[ n and the minimum temperature 
Tmin (the value below which the wind temperature is not allowed 
to drop). 

3.6. The adjustment of the wind formalism to ISA-Wind 

To be able to apply the analytical wind solution of our EOM dT6b 
to model a stellar wind from a given central star (with fixed 
stellar parameters) by using ISA- Wind to find numerically the 
unique solution, we had to adjust our wind formalism, i.e. our 
more accurate EOM, to the assumed EOM and wind velocity 
structure in ISA- Wind. The EOM ( [Tol l (and therefore our exact 
analytical wind solution) considers (allows) a line acceleration 
throughout the whole wind regime, starting above the radius r^ 6 , 
whereas the different EOM in ISA- Wind is only solved in the 
subsonic region by neglecting the line force. However then, ISA- 
Wind 'switches on' the line force somewhere below the connec- 
tion radius r con in the subsonic region by assuming a /3 velocity 
law above r con in the supersonic region. 

This inconsistency through the use of ISA- Wind (compared 
to our model assumptions and solutions) has been eliminated by 
introducing the parameter ro into our formalism for which the 
line radiation force is zero at radius fl' 6 '. Then, the final value of 
ro, together with the other remaining line acceleration or wind 
parameters, can be determined by fitting Eq. (l48l and the itera- 
tion procedure. 

Further, since ISA- Wind begins its computations already be- 
low (however close to) the stellar (i.e. photospheric) radius, all 
formulae derived in Sect.|2]have thus been applied with reference 
to the inner boundary (core) radius R m from where the numerical 
calculations of the wind model start. Therefore, the dimension- 
less variable of distance f here refers to R — R m (see the wind 
model for a typical O main-sequence star in Sect. [4] and the as- 
sociated diagrams in Sect. |2jl. 

3. 7. The chosen boundary values in the wind model 

In our following numerical wind model the inner boundary ra- 
dius has been chosen (constant throughout the whole iteration 
process) to be R m = 11 .757 Rq, situated at a fixed Rosseland op- 
tical depth of tr = 23. This corresponds then to a photospheric 
(stellar) radius of R p hot = 11 .828 Rq (defined as where the ther- 
mal optical depth is r e = 1 / y3) and an inner boundary density 
of pi„ = 1 .398 x 10~ 8 g/cm 3 at R in . 

This particular chosen fixed value of tr (or corresponding 
Pi„) at each step of the iteration cycle has no effect on the final 
wind parameters, i.e. in particular on the converged values of M 
and Voo, as further additional test iterations for a lower bound- 
ary value of T R < 23 (at the same /?;„) have shown. The reason 
why the exact value of p m is irrelevant, is, because we ensure 
by our iteration scheme that the sonic point, determined by the 
ISA- Wind computations, is at each iteration step n always lo- 
cated at the right radius, namely at the present critical radius of 
our EOM ( [TBT l and exact analytical solution. This current criti- 
cal radius is valid for all the present estimates of updated wind 
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parameter values obtained from the g'™ 6 curve fitting, i.e. in par- 
ticular for M„ and v Mn (see further explanations in Sects. 12.5.81 
andl!4l. 

3.8. General aspects of the boundary conditions for the 
uniqueness of the solutions 

To be able to solve the system of equations involving both the 
equation of continuity (0 and the momentum equation (fT0T > 
uniquely, one boundary condition in each case is required. To 
solve the momentum equation for the velocity v(r), the den- 
sity variable p was first eliminated by means of the equation 
of mass conservation, to obtain the density independent equa- 
tion of motion. This equation could then be solved uniquely for 
(e.g.) the wind velocity, because we demand a continuous and 
smooth trans-sonic flow through the critical point f c of Eq. ( [TBI 
for a stellar wind as one boundary condition. However, to be 
able to also solve the equation of continuity (0 uniquely for the 
mass-loss rate M by integration, as in Eq. ([8]), knowledge of the 
density p (r) for at least one local point r is required as another 
physical constraint. Thus, the density stratification for a given 
stellar wind (and therefore M) cannot be determined by merely 
using Eq. (O with the v (r) solution of the equation of motion. 

In more general terms, it can be shown that all solutions 
(Voo,M) lie on a curve, which fulfills the equation 

M v^ = const , with q w 2 , (54) 

since l/2Mv^ corresponds to the total wind energy (per time 
interval, i.e. equal to the kinetic wind energy far away from the 
potential of the central star). Hence, by providing an additional 
boundary condition, as, e.g., a given observed value for v„, the 
mass-loss rate would then be uniquely determined by a point on 
this curve (for which M = const v«?). 

However, since we are interested in self-consistent solutions, 
i.e. theoretical predictions of v(r) and M for a certain wind 
model from a given central star, without any prior knowledge 
of the terminal velocity v M , we follow the additional constraint 
that R (i.e. the inner boundary radius R[ n ) is situated at a fixed 
chosen optical depth t(R), or corresponding fixed density p(R) 
(cf. Sect. 13.7) . The latter value is given by Eq. (0 and the values 
for M and v (R), i.e. v; n . 

In our iteration process and by the use of ISA- Wind this cri- 
terion is always fulfilled at each iteration step. It is used to update 
(adapt) the value of v (R) at the inner boundary radius of the ve- 
locity assumed in ISA- Wind, by means of the new estimate of 
M obtained by MC-Wind (and Eq.[52]that depends on the value 
of Voo). After all wind parameters have converged, the final posi- 
tion of the sonic radius r s , i.e. the critical radius r c belonging to 
the point through which the unique final wind solution goes, is 
known. 

Hence, to summarise the mathematical point of view, the sec- 
ond necessary boundary condition we have used to be able to 
solve our system of equations (consisting of Eqs. l7land[T0l. or 
in other words, to find the unique values of M and Voo, is to pre- 
scribe the location of the critical point to be at the right position, 
i.e. at the radius r s of the sonic point. Whereas the density p (R) 
at the inner boundary radius is kept at a fixed value. 

4. Application: Results for an O-star 

In this section, we apply our theoretical results of Sect.[2]and the 
iterative procedure described in Sect. 13.31 to compute the stel- 
lar wind parameters for a typical 05-V main-sequence star. The 



Table 1. The different starting values in the iterations A - C for 
the wind from an 05-V main-sequence star." 



Iteration 


Vo<, [km/s] 


log M [MqI yr] 


A 


2020.0 


-5.5 


B 


6000.0 


-7.0 


C 


2020.0 


-7.0 



" For the stellar parameters, see Sect. [4] 



fixed stellar parameters (see e.g. Martins et al. 2005 and refer- 
ences therein) are: 

- = 40000 K 

- T = 0.214 

" ^core = 1 1.757 Rq 

- M = 40.0M© 

- log (L/Lq) = 5.5 

- Solar metallicity 

From which an isothermal sound speed of 18.16 km/s is ob- 
tained. 

The iteration method is tested for convergence by three dif- 
ferent iteration cycles A-C with three different starting values 
for Voo and M, see Table [TJ The reason for choosing these partic- 
ular values is to start the iteration process for Voo and M values 
that are either significantly higher, or lower, or lower for both, 
than the expected values. 

The results of each iteration step for each of the three itera- 
tion cycles are listed in Tables [2]to|4] In addition, some of these 
results are graphically displayed in Fig.[9]to[TT] where the val- 
ues of Voo, jS, r s , and log M, are plotted against the associated 
iteration step. 

From this, we can state that for each of the three runs all 
stellar and wind parameters converge to the same values (within 
statistical fluctuations) within approximately 10 iteration steps. 
However, to achieve a higher precision it is partly necessary to 
perform up to 15 steps in one iteration (see Tables [3] and |4] of 
iteration B and C). 

At the end of each iteration process A-C, we obtain the same 
converged radiative line acceleration g j.™ (r) (see def. Eq. [121 
of our stellar wind model as, e.g., in case of iteration A, plot- 
ted versus radial distance f in Fig. [TJ The dots represent the re- 
sults from the numerical calculation of g j.™ (r,) for discrete radial 
grid points r,. These numerical values are fitted well to the line 
acceleration function in Eq. ([T4l by the displayed fitting curve 
(solid line) with the (fitting, i.e. line) parameters go = 17661, 
y = 0.4758 ± 0.0002, 6 = 0.6878 ± 0.0003 and r = 1.0016 (ac- 
cording to Voo = 3232 km/s) at the end of iteration cycle A. to find 
the unique numerical solution The remaining converged stellar 
wind parameters (obtained by the same iteration cycle) as well 
as the resulting parameters from both other iterations B and C 
are listed in Table [5] 

Convergence is achieved when both values of the termi- 
nal velocity v M fl t and v*,, one obtained by the fitting formula 
Eq. d48l and the other by Eq. ( 1531 , have approached each other 
arbitrarily close (up to a desired precision). This is the case when 
the critical radius r c (determined by Eq.fTBl has become equal to 
the sonic radius f s (determined by ISA- Wind). 

Finally, we can summarise the numerical results of all iter- 
ations A-C by the mean values in Table [5] logM = -6.046 ± 
0.006 M©/yr, Voo = 3240 ± 37 km/s, f3 = 0.731 ± 0.005, or 
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Table 2. Iteration A for the wind from an 05-V main-sequence star. The variable stellar and wind parameters at each iteration step 
until convergence. 



Step Voo 


lOg M Voo, fil /3 


Tflt 


5m 


?0,flt 






no. [km/s] 


[M /yr] [km/s] 













-1 


2020 


-5.500 




1.0000 















6403 


-5.641 


2365 


0.8664 


0.7329 


0.4917 


1.0008 


1.0066 


1.0175 


1 


4022 


-6.194 


3846 


0.7981 


0.5963 


0.7048 


0.9999 


1.0087 


1.0101 


2 


3146 


-6.278 


4106 


0.7598 


0.5195 


0.8036 


1.0020 


1.0099 


1.0113 


3 


3003 


-6.226 


3909 


0.7355 


0.4711 


0.7747 


1.0030 


1.0101 


1.0115 


4 


2870 


-6.181 


3760 


0.7186 


0.4372 


0.7381 


1.0027 


1.0100 


1.0113 


5 


2827 


-6.141 


3684 


0.7047 


0.4094 


0.6882 


1.0024 


1.0100 


1.0111 


6 


2603 


-6.119 


3495 


0.7112 


0.4224 


0.7203 


1.0003 


1.0100 


1.0110 


7 


2450 


-6.078 


3301 


0.7352 


0.4704 


0.7349 


0.9940 


1.0099 


1.0111 


8 


2652 


-6.038 


3124 


0.7538 


0.5075 


0.7290 


0.9939 


1.0096 


1.0112 


9 


3211 


-6.047 


3232 


0.7379 


0.4758 


0.6878 


1.0016 


1.0095 


1.0112 



8 For the fixed stellar parameters L, T c a, R, M, and T, see Sect.[4] The line acceleration parameters yfi t (or equivalently /?), df, t and ro,fi t , and the 
terminal velocity Vco.sti were determined by the fitting formula Eq. d48t (or Eq.|49l>, applied to the results from a numerical calculation of the line 
acceleration gj.™ (?,). The parameter f' (in the /? velocity law, Eq. l41b and the sonic radius r s are output values from ISA-Wind, whereas log M is 
the improved estimated mass-loss rate numerically obtained (by MC-Wind and Eq.l52b. At each iteration step, the value of Va, was calculated by 
Eq. d53t and used as the new input value for the terminal velocity in the next iteration step, together with the new estimates of M and j3 (cf. 
description of iteration process in Sect. l3.3t . Convergence is achieved when the values of Voo.at and v„ have become equal (i.e. here, close) to each 
other. In this case, the condition that the critical radius r c (determined by Eq.l 1 8t has to be equal to the sonic radius r s = 1.01 1 (determined by 
ISA- Wind) is fulfilled. 



Table 3. Iteration B for the wind from an 05-V main-sequence star. The variable stellar and wind parameters at each iteration step 
until convergence. 



Step 


Voo 


log M 


Voo, fit 


P 


Tat 


5fit 


fo,Rt 






no. 


[km/s] 


[M Q /yr] 


[km/s] 














-1 


6000 


-7.000 




1.0000 















4208 


-6.863 


7011 


1.0058 


1.0116 


1.3468 


0.9976 


1.0107 


1.0142 


1 


4396 


-6.655 


5511 


0.9347 


0.8695 


1.1506 


1.0042 


1.0100 


1.0153 


2 


3683 


-6.561 


4968 


0.9020 


0.8040 


1.1385 


1.0002 


1.0101 


1.0135 


3 


3861 


-6.436 


4588 


0.8275 


0.6550 


0.9511 


1.0051 


1.0099 


1.0133 


4 


3291 


-6.385 


4432 


0.7843 


0.5687 


0.9136 


1.0039 


1.0102 


1.0121 


5 


3116 


-6.297 


4152 


0.7535 


0.5070 


0.8314 


1.0036 


1.0102 


1.0118 


6 


2849 


-6.230 


3922 


0.7416 


0.4833 


0.7919 


1.0014 


1.0101 


1.0115 


7 


2716 


-6.165 


3596 


0.7422 


0.4844 


0.7875 


0.9996 


1.0100 


1.0114 


8 


3311 


-6.117 


3544 


0.7118 


0.4237 


0.6827 


1.0046 


1.0098 


1.0113 


9 


2768 


-6.163 


3765 


0.6948 


0.3897 


0.6935 


1.0031 


1.0100 


1.0109 


10 


2321 


-6.120 


3460 


0.7213 


0.4426 


0.7561 


0.9953 


1.0101 


1.0110 


11 


2922 


-6.041 


3228 


0.7198 


0.4396 


0.6698 


1.0008 


1.0097 


1.0112 


12 


2629 


-6.080 


3385 


0.7286 


0.4572 


0.7009 


0.9968 


1.0097 


1.0108 


13 


3022 


-6.063 


3358 


0.7168 


0.4336 


0.6633 


1.0018 


1.0097 


1.0111 


14 


2493 


-6.102 


3397 


0.7334 


0.4669 


0.7447 


0.9955 


1.0098 


1.0108 


15 


3008 


-6.056 


3307 


0.7225 


0.4449 


0.6696 


1.0013 


1.0097 


1.0112 



8 See description of parameters, analogous to iteration A, in Table[2]above. 



y = 0.462±0.009, 5 = 0.6811+0.0001 and h = 1.0014±0.0001 
with an average critical (i.e. sonic) radius of r c = 1.0098. 

The model results for the wind velocity v (r) as a function 
of radial distance r from the typical 05-V-star, with the above 
mentioned converged wind and line acceleration parameters of 
iteration A, are displayed in Figs. [5]- [7] and have been already 
discussed in Sect. 12.51 

Figure [5] shows the comparison of the approximated and the 
exact wind solution with the converged line acceleration param- 
eters go, y, 5 and fo of iteration A in Eq. (1141) . where the approxi- 
mated wind curve was obtained using Eq. (1391 , whereas the exact 
wind solution was obtained from Eqs. ( [30) , ( f3TT >, and d33l ), with 



a critical radius of r c = r s — 1.0110. A more detailed compar- 
ison of the approximated and the exact wind solution is shown 
in Fig.|7j where both curves are displayed for three different ra- 
dial ranges. The only noticeable discrepancy between the wind 
velocity curves is in the subsonic and sonic region where the ap- 
proximated wind solution, derived for (and only supposed to be 
applied in) the supersonic region, is no longer valid. Figure [6] 
displays the same approximated wind solution together with 
the f3 wind velocity curve, where the latter was obtained from 
Eq. ( f4TT > with the converged values for the exponent /3 = 0.7379 
and the terminal velocity Voo = 3232 km/s of iteration A (with 



14 Patrick E. Miiller and Jorick S. Vink: Radiation-driven winds from massive stars 

Table 4. Iteration C for the wind from an 05-V main-sequence star: The variable stellar and wind parameters at each iteration step 
until convergence. 



Step Voo 


log M Voo, « /3 


Tflt 


5m 








no. [km/s] 


[M©/yr] [km/s] 













-1 


2020 


-7.000 




1 .0000 















4161 


-6.448 


4467 


0.9230 


0.8459 


1.0278 


1.0052 


1.0097 


1.0203 


1 


4080 


-6.430 


4506 


0.8463 


0.6927 


0.9593 


1 .0043 


1.0097 


1.0130 


2 


3431 


-6.407 


4515 


r\ Ton 

0.7937 


0.5875 


0.9191 


1 .0040 


1.0101 


1.0121 


3 


3340 


-6.326 


4260 


0.7553 


0.5107 


0.8416 


1.0052 


1.0102 


1.0119 


4 


2791 


-6.272 


4087 


0.7300 


0.4600 


0.7965 


1.0025 


1.0102 


1.0115 


5 


2619 


-6.179 


3683 


0.7321 


0.4642 


0.7760 


0.9995 


1.0101 


1.0114 


6 


2924 


-6.111 


3484 


0.7251 


0.4501 


0.7078 


1.0015 


1.0099 


1.0113 


7 


2917 


-6.116 


3509 


0.7098 


0.4196 


0.6994 


1.0028 


1.0098 


1.0110 


8 


2752 


-6.116 


3560 


0.7073 


0.4146 


0.6899 


1.0014 


1.0099 


1.0110 


9 


2688 


-6.094 


3450 


0.7084 


0.4169 


0.6930 


1.0007 


1.0099 


1.0110 


10 


2800 


-6.076 


3375 


0.7052 


0.4105 


0.6833 


1.0019 


1.0098 


1.0110 


1 1 


2828 


-6.082 


3438 


0.6978 


0.3956 


0.6657 


1.0024 


1.0098 


1.0109 


12 


2416 


-6.090 


3369 


0.7119 


0.4238 


0.7133 


0.9969 


1.0099 


1.0108 


13 


2572 


-6.039 


3142 


0.7394 


0.4788 


0.7235 


0.9950 


1.0097 


1.0111 


14 


3087 


-6.035 


3181 


0.7320 


0.4639 


0.6859 


1.0012 


1.0095 


1.0111 



" See description of parameters, analogous to iteration A, in Table|2]above. 



Table 5. The results of the different iterations A-C for the wind from an 05-V main-sequence star. The converged values of the 
wind parameters v,*,, ft, y, 5 ,fo in the radiative line acceleration Eq. d48l i (or Eq.|49] respectively), and mass-loss rate log M of each 
iteration cycle." 



Iteration 


log M 




P 


7 


6 


Hi 


h 




[M /yr] 


[km/s] 












A 


-6.047 


3232 


0.7379 


0.4758 


0.6878 


1.0016 


1.0110 


B 


-6.056 


3307 


0.7225 


0.4449 


0.6696 


1.0013 


1.0083 


C 


-6.035 


3181 


0.7320 


0.4639 


0.6859 


1.0012 


1.0101 


Mean values 


-6.046 ± 0.006 


3240 ± 37 


0.731 ± 0.005 


0.462 ± 0.009 


0.6811 ± 0.0001 


1.0014 + 0.0001 


1.0098 



" In addition, the critical radius r c (equal to the sonic radius r s ), determined by Eq. J 1 8b . is displayed in the last column, as it is needed as a further 
parameter to evaluate the wind solution by Eqs. J30t and d31t . The mean parameter values (with error estimations) have been calculated from the 
results of all iterations A-C. 



K = 1.0095 as used in IS A- Wind). Both curves approach each 
other and the same velocity limit at large radial distances r. 

These comparisons highlight the good agreement between 
all velocity curves in the supersonic region, except for the notice- 
able difference of the and approximated wind velocity curve 
at intermediate distances r > 5 from the star (see Fig. [6]). In 
the supersonic region, the /3 velocity curve always lies (e.g. at 
radial distance of r — 10 by 12%, and at f — 100 only by 3%) 
above the approximated solution. Compared to the exact wind 
solution, this means that the /?-wind law overestimates the real 
wind velocity over the entire range of radial distances. 

This difference at intermediate distances is mainly due to our 
approximation (Eq. |4"7T i for the relationship between the expo- 
nents p and y (cf. the comparison of the /3 law with our approx- 
imated velocity) that is only approximately independent of the 
exponent 5, if the latter is close to 1 . Since in our wind model 
the value of 6 has actually converged to 0.7, this approximated 
relation does not provide a better value for f3 to adapt the /3 curve 
closer to the approximated (and also real) velocity curve in the 
middle wind regime. Nevertheless, we first had to introduce the 
additional exponent 5 in our line acceleration function, Eq. ( fT4l . 



to gain an optimal fitting to the numerical calculation of gj.™ (r) 
at those middle radial ranges of the wind. 



5. Discussion and conclusions 

In this paper, we have derived analytical solutions for the ve- 
locity structure for any mass outflow or inflow situation, through 
use of the Lambert W-function. For the case of a radiation-driven 
wind, we described the line acceleration as a function of stellar 
radius g(r), and therefore found, as a mathematical consequence, 
the critical point of our equation of motion to be the sonic point 
(in contrast to the critical point of the different EOM of CAK). 

The reason why the line acceleration can be expressed as 
a function of radius only can be justified as follows: By solving 
the basic hydrodynamic equations under the assumption of a sta- 
tionary (i.e. steady) spherical flow, we assume that the problem 
is solvable time-independently and all physical quantities, in par- 
ticular the solution functions, and especially the line acceleration 
must be expressable as a function of radius only. This occurs at 
the latest after the physical process has reached an equilibrium. 
The evidence is given by this work, where, at the end of three 
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Fig. 9. Iteration A for the wind from an 05-V main-sequence star: The varying wind parameters Voo.fit ( see upper left diagram) 
and /3 (upper right), both obtained by the fitting to numerical results for the line acceleration, the numerically determined sonic 
radius r s (by IS A- Wind, see lower left), and the stellar mass-loss rate log M (calculated by MC-Wind, see lower right), plotted 
versus the step numbers of iteration A until convergence. For the fixed stellar parameters L, T e ff, R, M and F, see Sect. [4] See 
also the description of parameters in Table [2] 



different iteration cycles, g(r) converges to the same function of 
r for the same stellar wind model. 

The numerical calculation of g(r) was performed through 
Monte Carlo simulations and the wind parameters were simul- 
taneously solved in an iterative way for a typical 05-V star. 
Our computations converged to logM = -6.05 MQ/yr, v M = 
3240 km/s, /? = 0.73 and r s = 1.011 (or expressed equivalently 
to go = 17392, y = 0.46 and 6 = 0.68). 

The question is how these self-consistently derived wind pa- 
rameters compare to earlier computations and observational con- 
straints. First, when we employ the mass-loss recipe of Vink 
et al. (2000) using the stellar parameters of r eff = 40000 K, 
M = 40.0 Mq, and log (L/Lq) = 5.5, for an unspecified termi- 
nal wind velocity, so that the recipe uses the standard value of 
Voa/v esc of 2.6, the recipe yields a value for the mass-loss rate of 
logM = -5.9 MQ/yr. In other words, our self-consistent com- 
putation yields a mass-loss rate in agreement with the Vink et al. 
recipe. 

We now discuss the values that we found for the velocity 
field. First of all there is the shape of the velocity field, y = 0.46 
and 6 = 0.68, or more widely recognised as /3 — 0.73. This 
value still agrees with earlier modified CAK calculations (e.g. 
Pauldrach et al. 1 19861 1 as well as empirical line profile analyses 
(e.g. Puis et al. [1996l Mokiem et al. 1200751 ). 

We finally turn our attention to our derived terminal wind 
velocity of v M = 3240 km/s. This value is still in line with ob- 
servations, where blue edges of P Cygni profiles do not extend 
beyond ~ 3000 km/s for O star dwarfs. We assumed the same 



solar abundances as in Vink et al. ( 120001 1, which are based on 
those of Allen ( I1973I ). However, a recent 3D analysis of the so- 
lar atmosphere by Asplund ( 2005 ) suggests that the solar carbon, 
nitrogen, and oxygen abundances may be lower than we have as- 
sumed. As these intermediate-mass elements are important line 
drivers in the supersonic portion of the stellar wind (Vink et al. 
1999, Puis et al. 2000), this may have implications for the outer 
wind acceleration, and hence might somewhat reduce the termi- 
nal wind velocity we obtain. 

Our current models are non-rotating, while real O stars rotate 
rather rapidly. Friend & Abbott (1986) showed that the termi- 
nal wind velocity drops with equatorial rotation speed according 
to Voo/vesc K (1 !iot_)0.35_ ^ s q stars ro t a t; e a j a significant 

^'break-up 

fraction of their break-up speed, our non-rotating models pro- 
vide just the upper limit to the terminal wind velocity (and the 
lower limit to the mass-loss rate). Further investigations, taking 
into account the effects of rapidly rotating O stars, predict (up to 
~10%) lower terminal velocities. 

The above-mentioned issues will form the basis for future 
investigations, for which we have laid the groundwork in the 
present paper. 
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Fig. 10. Iteration B for the wind from an 05-V main-sequence star: The varying stellar and wind parameters Voo.at ( see upper left 
diagram),/? (upper right), both obtained by the fitting to numerical results for the line acceleration, the sonic radius r s (determined 
by ISA-Wind, see lower left), and the mass-loss rate log M (calculated by MC-Wind, see lower right), plotted versus the step 
numbers of iteration B until convergence. 
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Fig. 11. Iteration C for the wind from an 05-V main-sequence star: The varying stellar and wind parameters Voo.a ( see upper left 
diagram),/? (upper right), both obtained by the fitting to numerical results for the line acceleration, the sonic radius r s (determined 
by ISA-Wind, see lower left), and the mass-loss rate log M (calculated by MC-Wind, see lower right), plotted versus the step 
numbers of iteration C until convergence. 



